Stability of strained heteroepitaxial systems in (1 + 1) dimensions 
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We present a simple analytical model for the determination of the stable phases of strained 
heteroepitaxial systems in (1 + 1) dimensions. In order for this model to be consistent with a 
subsequent dynamic treatment, all expressions are adjusted to an atomistic Lennard- Jones system. 
Good agreement is obtained when the total energy is assumed to consist of two contributions: the 
surface energy and the elastic energy. As a result, we determine the stable phases as a function 
of the main "control parameters" (binding energies, coverage and lattice mismatch). We find that 
there exists no set of parameters leading to an array of islands as a stable configuration. We however 
show that a slight modification of the model can lead to the formation of stable arrays of islands. 

PACS numbers: 68.35.-p, 68.43.Hn, 68.65.-k, 68.65.Hb 



I. INTRODUCTION 

It appears today that self-assembly is not only one of 
the most elegant avenues for the production of devices 
based on quantum dots (QD), but also one of the most 
promising. Basic understanding of the physics of the for- 
mation of arrays of islands should ultimately lead to the 
realization of such exciting concepts as spintronicsi*^ and 
quantum dot cellular automata, 3 Driven by these possible 
developments, considerable effort has been devoted to un- 
derstanding and predicting the conditions necessary for 
ensuring the stability of arrays of islands. Simple argu- 
ments based on the scaling of the energy as a function of 
the volume of the islands^ indicate that any system natu- 
rally undergoes ripening and, therefore, the only relevant 
observation is the (very long) time scale needed for the 
system to reach equilibrium^ It is becoming clear, how- 
ever, that a realistic energy function can lead to arrays of 
islands as equilibrium configur at ions ; 4 i 6 1 7 1 8 1 9 as can also 
be deduced from experiment (see, e.g., Ref. IT(il for a re- 
view) . 

The dynamics of formation of arrays of islands 
has to some extent been investigated using atomistic 
modelSfii*i2iiii4 but is not yet fully understood. For 
instance, the importance of nucleation in the early stage 
of array formation is still unclear^ The long-term goal 
of our work is to address such questions and to provide 
a coherent picture of island formation in heteroepitaxial 
systems, duly taking into account changes in the energy 
landscape arising from the lattice mismatch between the 
two components of the system. We aim to achieve this 
using a kinetic Monte Carlo (KMC) model, whereby the 
particles evolve according to the relative probabilities for 
hopping from one site to a neighbouring site on a fixed 
lattice. The main difficulty of this approach, which we 
are still in the process of developing, resides in properly 
modulating the energy barriers to account for elastic con- 



tributions generated by the lattice mismatch between the 
two species. We note that this problem could in princi- 
ple be approached using molecular dynamics (MD). How- 
ever, the timescales involved in the present problem are 
such that MD calculations are not feasible at this time. 

For lack of an accurate model for the kinetics of is- 
land formation, there is definite interest in the investiga- 
tion of the equilibrium properties of heteroepitaxial sys- 
tems, i.e., the surface morphology as determined by the 
various parameters that control the physics of the sys- 
tems (lattice mismatch, coverage, binding energies, etc.). 
Our objective here is to present and discuss a continu- 
ous model suited for this purpose. This will enable us, 
in particular, to identify the regions of parameter space 
where the formation of stable, coherent arrays of island 
are favored. We develop an analytical expression for the 
zero-temperature total energy of a strained array of is- 
lands. To ensure consistency, an important constraint 
on this model is that it should be expressed in terms of 
quantities that can be "exported" to a subsequent KMC 
model which will allow the dynamics to be investigated. 
We do this by assuming Lennard- Jones interactions be- 
tween the two types of atoms and considering a triangular 
"(1 + l)-dimensional" geometry. This defines a reference 
system (which we will call the L J system) on which both 
the static and the dynamic models should be mapped. 
Note that we use the term "(1 + l)-dimensional" , rather 
than "2-dimensional" , to indicate that one of the spatial 
dimensions is height (the z component), not to be con- 
fused with the usual two-dimensional case where atoms 
move in the xy plane. 

We thus obtain an expression for the difference in en- 
ergy AE between a system with a flat layer of adsorbed 
atoms and a system where islands have formed. For a 
given set of control parameters, this quantity is mini- 
mized with respect to the size of the islands, the dis- 
tance between the islands, and the thickness of the wet- 
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ting layer, so as to determine the equilibrium state of 
the system. We find that good agreement between the 
LJ system and the continuous model can be achieved by 
considering only two contributions in AE: the surface en- 
ergy and the elastic energy, the latter arising from both 
relaxation and substrate-mediated island-island interac- 
tions. This general approach is in many respects similar 
to that proposed by Combe et al. 4 (CJB). Apart from 
the fact that our model allows the lattice misfit and the 
energy parameters to vary, the main differences lie in the 
details of the method, as discussed below. 

An important conclusion of our work is that no single 
set of control parameters leads to an array of islands as 
a stable configuration. This is in contradiction with the 
results of CJBji the discrepancy arises from the arbitrari- 
ness in the choice of the parameter (z ) which represents 
the characteristic length for the decay of the adsorption 
energy above the surface. If we relax the constraint on 
the choice of this parameter (i.e., consistency with the LJ 
model is not imposed), then stable arrays of islands are 
found, along with a "cracks" phase, made of flat islands 
that touch at their base. 



In practical calculations, the LJ interaction must be 
cutoff at some distance r c . The choice of r c affects slightly 
the physical properties of the system, notably the equi- 
librium interatomic distance r cq , the cohesive energy uq, 
and the elastic constants. (See Section III B 21 for de- 
tails on the calculation of the elastic constants). Table 
[3 presents the dependence of these important quantities 
on the cutoff radius. 



r c 


r eq (<x) 


uo (e) 


H = A (e) 


1.000 (1) 


1.1225 


-3.000 


31.18 


1.732 (2) 


1.1159 


-3.222 


33.48 


2.000 (3) 


1.1132 


-3.319 


34.49 


2.646 (4) 


1.1122 


-3.356 


34.87 


3.000 (5) 


1.1119 


-3.364 


34.96 


oo 


1.1115 


-3.382 


35.15 



TABLE I: Computed values of some important quantities as 
a function of the cutoff distance r c (and number of nearest- 
neighbor shells): r cq is the equilibrium interatomic spacing 
in units of a, ito is the cohesive energy in units of e, and fj, 
and A are the two Lame parameters (both are equal in this 
geometry). 



II. THE MODEL 



A. The LJ system 



The Continuous Model 



As mentioned above, the reference system consists of 
atoms occupying the sites of a (1 + l)-dimensional trian- 
gular lattice and interacting via the Lennard-Jones po- 
tential: 



U(r) = 4e 



r ■ 



r 



(1) 



The atoms are of two types: substrate (S) and adsorbate 
(A). Thus, there are three different types of interactions 
and six different Lennard-Jones parameters that define 
the total energy: {ass, &aa, o~sa} and {ess, £aa, ssa}- 
Since one of the cr's and one of the e's fix the length scale 
and energy scale, respectively, there are only four free 
parameters. This number can be reduced to three by 
applying the Lorentz-Bcrthelot combination rule: 



as A 



1/ 



o~aa) 



(2) 



Hence, there is a single degree of freedom as far as length 
scales are concerned, which can be expressed in terms of 
the lattice mismatch a, defined as 



a A A - ass 
&ss 



(3) 



The other two degrees of freedom are the binding ener- 
gies between adsorbate atoms, £aa, and between adsor- 
bate and substrate atoms, esA- These parameters are 
independant and 65^4 > 6aa is the wetting condition. 



As discussed in the Introduction, our purpose is to 
evaluate the energy difference between a system in which 
the adsorbate atoms form islands and one in which they 
form a uniform layer on top of the substrate: 



AE 



E; 



layer • 



(4) 



This energy difference can be decomposed into surface 
and elastic contributions, 

AE = AE SUT f acc + AS c i as tic , (5) 

and is a function of the following parameters: 

• a, the lattice mismatch; 

• ess, the binding energy between two atoms of the 
substrate; 

• caa, the binding energy between two atoms of the 
adsorbate; 

• esA, the binding energy between an atom of the 
substrate and an atom of the adsorbate; 

• 8, the coverage, expressed in monolayers (ML); 

• h, the height of the islands, expressed in ML; 

• L, the width of the islands at their base; 

• z, the thickness of the wetting layer, expressed in 
ML; 
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FIG. 1: Geometry of the system (a) with islands and (b) 
without islands. The shaded region is the substrate and the 
white region is the adsorbate. Islands are assumed to have 
the shape of an isosceles trapezoid, with contact angle j. 



FIG. 2: Displacements of the atoms of a relaxed island (L = 
20, h — 10). (a) tensile (a — —1%); (b) compressive (a = 
+ 1%). The lengths of the arrows are about 30 times the 
actual atomic displacements. 



• d, the distance between the centers of two neigh- 
boring islands. 

As mentioned earlier, one of the binding energies, say 
esS) fixes the energy scale. The last five parameters de- 
scribe the geometry of the system (cf. Fig. QJ; they are 
integers but we will assume that they are real in order to 
facilitate the calculations. Since h, 6 and z are expressed 
in ML, the actual height is obtained by multiplying by 
the thickness of one monolayer, ^r eq . 

The conservation of atoms (or volume) between the 
two configurations imposes a constraint on the geometric 
parameters: 



9d = zd + h{L-h/2). 



(6) 



Since a, ess, csa, £AA and 9 are assumed to be known 
from experiment for a given material, i.e., they can be 
considered as control parameters, the energy need only 
be minimised with respect to the three parameters L, h 
and z, d being determined by the constraint ifrjfl. 

In the next two subsections we develop expressions for 
the two energy contributions in Eq. JSJ as a function of 
the various parameters of the problem. They are derived 
largely from theoretical considerations but, in some cases, 
ad hoc terms are introduced in order to ensure that the 
model is consistent with the LJ system; in these cases, 
we proceed as follows: first, we generate a set of con- 
figurations with some number of adsorbate atoms placed 
in the shape of a trapezoid. Then, with suitably chosen 
parameters (e and a), we let the whole system (includ- 
ing a substrate large enough to make finite size effects 
negligible) relax to minimum energy. Periodic boundary 
conditions are used along the x direction. Figure |21 illus- 
trates the atomic displacement of the atoms of an island 
as a result of relaxation. 



1. The Surface Energy 

We first determine the adsorption energy of an island 
of type S (i.e., "substrate on substrate") whose actual 
height is h(x), or h(x) when expressed in ML (see Fig. [3}- 
This energy is the sum of bulk and surface contributions: 



^uoess + {L — L)jss, 



(7) 



where V is the volume of the island, it ess the cohesive 
energy per atom, [L — L) the increase of surface due to 
the island, and 7ss the surface energy density. Using L = 
J Vdx 2 + dz 2 = J y/T+W(x)dx, Ef and can readily be 
rewritten in terms of h(x): 



hu ess 




1-1 7ss 



dx, (8) 



where the factor | comes from the substitution of h(x) 
by &h(x). 
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FIG. 3: Island of shape given by h(a 
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The surface energy density jss is proportional to ess', 
we can write 



7SS = Cess 
Iaa — Ceaa 

Isa — C(ess + £aa — 2£sa), 



(9) 



where C is a constant that depends on the number of 
neighbors taken into account in the model. 

If the island is of type A (i.e., adsorbate on substrate), 
now, the adsorption energy is easily obtained from Eq. 
101 by substituting caa and jaa for ess and 7ss and 
adding a term describing the interaction between the is- 
land and the substrate. If we assume nearest-neighbor 
interactions, this term involves only the atoms at the in- 
terface between the island and the substrate. In a more 
general situation, the adsorption energy for an atom at 
position z above the surface (in ML) can be written 



E ad (z) 



^ad + ( e AA ~ € S a)v(z), 



(10) 



where E® d is the "generic" adsorption energy for the case 
6sa = £aa and rj(z) is a function which decreases with 
z, with characteristic length zq. As in similar works^ii& 
we assume the form: 



r)(z) = Ae- z ' zo - 



(11) 



the parameters A and zq can be determined by fitting 
to the total energy of a particular atomic model (here 
Lennard- Jones). Using (|1U[I. we find that the adsorption 
energy of a vertical column of h atoms is proportional to 



3=1 

Hence, we obtain: 



island 




I _ g-?i/ 2 o 
e l/z _ 1 



\h> 2 + 1 - 1 J Iaa 



(12) 



hu €AA 

+B{e A A-esA)(l-e- h / z °) dx, (13) 



E 



where B = Aj (e 1 ' z ° — 1). It can easily be shown that B = 
2C by taking the limiting case of a very thick adsorbate. 
Overall, therefore, only two parameters need to be fitted 
to the atomic model, namely B and zq. This fit was 
carried out on systems with uniform adsorbed layers of 
thickness 9 [i.e., h(x) = 6 => h'(x) = 0]. In this case, the 
energy difference between a system with an adsorbate of 
type A and an adsorbate of type S is 



Assuming h{x) describes a trapezoidal shape, we can 
finally write the first term of as 

A£ surfaco = 2C{e AA - e SA )[d - (d - L + z )e- z ' z ° 

-(L-h- zti)e-( h+z V zt) - d(l - e- e/z0 )} + Ce AA h. 

(15) 

2. Elasticity 

It is an interesting fact that the theory of elasticity, 
which deals with continuous media, can accurately de- 
scribe systems as small as a few tens or hundreds of atoms 
(see Ref. for instance). In this work, we exploit this 
property to construct (at least in part) the analytical 
expressions entering the second term of (JSJl. Unfortu- 
nately, we know of no analytical solution to the elastic- 
ity differential equations for a system with the boundary 
conditions illustrated in Fig. [IJ,). This difficulty can be 
circumvented by making some assumptions on the force 
distribution caused by the island on the substrate. As 
argued by Tersoff and TrompflS. when the island is very 
flat (h <C L), one can assume that there is no strain 
relaxation in the z direction. This leads to a force dis- 
tribution which is proportional to the derivative of the 
height function, / oc h'(x). In more general cases, how- 
ever, we found that a (truncated) linear force density was 
a better assumption. Since higher islands deform the sub- 
strate more efficiently (for a given volume), we used the 
following expression: 



kx if \x\ < -§ 
otherwise, 



(16) 



where k is a constant (to be determined) which depends 
on the lattice mismatch, the shape of the island and the 
elastic modulii of both the island and the substrate. 

Following the guidelines presented in Ref. ITfll we 
find that the Green's tensor for a semi-infinite two- 
dimensional plane is 



Gxx 
G X z 
G zx 

G zz 



(A + /x) x 2 - (A + 2/i) r 2 log (r) 

2-ir/j, (A + n) r 2 
(A + fi) xz — /ir 2 arctan (|) 

27T/i (A + (i) r 2 
(A + fi) xz + \ir 2 arctan (|) 

27T/x (A + [i) r 2 
- (A + n) x 2 - (A + 2/x) r 2 log (r) 
271-^4 (A + fi) r 2 : 



(17) 



E 



layer 



E l g ycl = d[8uo(e A A - ess) 



+ B(e A A-e S A)(l-e~ e/zo )}. (14) 

The numerical calculations yield B = 2.53 and zq = 0.39 
(see Appendix lAl for details). 



where r 2 = x 2 +z 2 . The x component of the displacement 



field on the surface of the substrate is found to be 



u x (x,z = 0) = j^Gxxix, z = 0)f(x - x')dx' 



(18) 
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where is a scale-independant function: 



2^+ (i-e log 



(19) 



with k = i6^ s anf l Ms i s one °f the Lame parameters of 
the substrate. 

Because the atoms lie on a triangular lattice and in- 
teract via a radial potential, the two Lame parameters fj, 
and A must be equal, given by 



(20) 



the sum being carried out over all lattice points within 
a given cutoff radius. The Young modulus and Poisson 
ratio for a 2-dimensional system arc 24 



E 



4/i(A + n ) 8 
A + 2/i 
A _ 1 
A + 2/x _ 3 



(21) 



With these assumptions, we can construct expressions 
for the two principal elastic contributions to the total 
energy arising from the presence of a mismatched island 
on a substrate, viz. the energy due to the mutual strain 
between the island and the substrate, and the island- 
island interaction energy mediated by the substrate: 



AE. 



clastic 



= AE stra 



E 



interaction- 



(22) 



The elastic energy per unit volume (more precisely unit 
area in the present case) of a strained uniform adsorbed 
layer is 



jglaycr 
clastic 

dV 



(23) 



Ea being the Young modulus of the adsorbate. Hence, a 
portion of width d of this ^#-thick layer has an energy 



E 



layer 
clastic 



Similarly, we can write 



E 



sland 
clastic 



x/3 



-E A a 2 v(L. h), 



(24) 



(25) 



where v(L,h) is a function (to be determined) having 
units of volume. CJE 4 write this as 



v(L,h) = R(r)V, 



(26) 



where r = h/L is the aspect ratio, V is the volume of the 
island, and R(r) is a dimensionless function bounded be- 
tween and 1. Evidently, for scaling reasons, any v(L,h) 
can be written in this form. In the present work, how- 
ever, we choose to determine v{L, h) without invoking 



this scaling ansatz. This choice is mainly dictated by ac- 
curacy considerations: in order to fit R{r) to numerical 
data, it is necessary to divide the computed elastic energy 
by the volume, thereby reducing the relative importance 
of large islands (see appendix iBl for more details). 

We find that an excellent fit to the Lennard-Jones data 
is obtained with 



v(L,h) 



y/3L 



2 r 



-ch/(L-h) 



(27) 



where c is the only (unitless) parameter to be adjusted 
to the data. Numerical calculations on the LJ system 
yield c = 13.5 (see Appendix lAl for details). Ratsch and 
Zangwill 20 have developed a similar function from theo- 
retical considerations; the only difference is the denomi- 
nator of the exponent, where they use L instead of L — h 
and find c = 2\/37r « 10.9. 

In summary we have, for the elastic energy due to the 
strain between the island and the substrate: 



AE S . 



V3 



fj, A a 



L 2 

c 



-ch/(L-h) 



zd — 9d 



(28) 

where wc include the contribution arising from the pres- 
ence of a wetting layer of thickness z in the system with 
islands. 

The substrate-mediated interaction energy between 
two islands a distance d apart — the second term in Eq. 
(|22|l — can be written as a surface integral, 



E, 



l-isl 



(d) 



u x {x)f(x - d)dx, 



(29) 



where u x (x) is the displacement field of the first island 
and / is the surface force distribution of the second is- 
land. Using 11 6|l and (|19|l . we get: 



Ek, 



l-isl 



(d) 



3L+ + 2d z L z +4dL d log 
d 2 -L 2 



pL 2 

+2(d 4 - 3d 2 L 2 ) log 



d-L 



d + L 
(30) 



Note that this expression remains finite when d — ► L. 
Since we are interested in the interaction energy for an 
array of islands, we have to sum up all contributions com- 
ing from the islands at position ...,—2d,—d,d,2d,... In 
order to get a simple expression for our model, we may 
replace i?isi-isi by the first few terms of its asymptotic 
series, 



Eis 



isl-isl 



2ir 2 [i s L 2 K 2 ( 2L 2 L 



3L b 



3nd 2 5nd 4 35nd 6 



The sum of all possible contributions can now be carried 
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out exactly: 



E, 



interaction 



oo 

E 

j=—oo 



E-, 



isl-isl 



(jd) 



47r 3 Ms £ 2 K 2 
81 



d 2 



n 2 L 4 
50d 4 



7T 4 L 6 \ 

1225d 6 J 



(32) 



We still have to determine the dependence of k [cf. Eq. 
(|19f) ] on the parameters of our model. Let us first recall 
that this quantity is proportional to the amplitude of the 
force distribution exerted on the substrate by the island 
and, therefore, it also determines the amplitude of the 
displacement field in the substrate. We assume k to have 
the following form: 



a k^ a ,Hs) K gco (h,L), 



(33) 



where k m depends only on the elastic coefficients of the 
susbtrate and the adsorbate, and K seo depends only on 
the geometry of the island, n is linear in a because Eq. 
(|19|l has itself been derived within the framework of the 
linear theory of elasticity. 

For we propose the ad hoc expression 



HA 



HA + HS 

which satisfies the three important limiting cases 



(34) 



HS > HA 
MS < HA 

MS = HA 



- 

• 1 
1 

2' 



(35) 



The first relation corresponds to the case of a substrate 
which is much more rigid than the adsorbate, the second 
is the opposite case, and the third is the case of equally 
rigid substrate and adsorbate. The geometry factor K geo 
can be determined by fitting to the LJ data. We find the 
following expression to yield good results: 



-aiL+tt2 



-b 1 h/(L-h)+b 2 



(36) 



The first factor, ^, is such that d x u x , the x-component 
of the substrate's strain tensor, is always between and 
a. The second factor ensures that K geo goes rapidly to 1 
when L is larger than a few atoms. The last factor is a 
scale invariant term which is maximum when the aspect 
ratio 4 is 1. After fitting the displacement field to the 
L J data (see Appendix EJ , we find the following set of 
values for the parameters entering the above expression: 



oi - 12.6/r eq , a 2 = 0.028, 
h = 0.033, b 2 = -1.35, 



(37) 



where r oq is the equilibrium interatomic distance (see pre- 
vious section). 

Putting everything together, we obtain the following 
long expression for the interaction energy in Eq. I|22(l : 



interaction 



47rVs£ 2 
81 



Ha 



4 [l A + fJbs 



1 -e 



— a\L-\-a2 



)('-- "™)] 2 (S + S? + ^)- < 38) 



We are now in position to construct the phase diagram of 
our model as a function of the various control parameters. 



III. RESULTS AND DISCUSSION 

A. Phase Diagrams 

The energy difference AE between a system in which 
the adsorbate atoms form islands and a system in which 
they form a uniform layer on top of the substrate, Eq. 
(@J, can be expressed as 



AE = A£L 



rfacc 



AE stIB 



E; 



interaction j 



(39) 



where the three terms are given by Eqs. Q15|). (|28[l. and 
(|38|l . respectively. We now proceed to determine the 
equilibrium configurations of the system as a function of 
the control parameters, viz. binding energies, mismatch 
and coverage. The parameter space is sampled as follows: 

e AA = 0.7, 0.8,..., 1.2, 1.3 

e SA = 0.7, 0.8,..., 1.2, 1.3 

a = 0%,1%,...,9%,10% 

6 = 1,...,14,15, 

and ess — 1 sets the energy scale. We took positive 
values of a only since AE depends quadratically on this 
parameter. For every possible combination of the above 
parameters, we determine h, L, d and z such that AE 
is minimal. If the minimum AE is larger than zero, the 
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0% 2% 4% 6% 8% 10% 

Lattice mismatch (a) 



FIG. 4: (Color online) A typical phase diagram; here eaa = 1 
and ssa = 1.1. The three phases are: ^ Frank- van der 
Merwe; ripening islands with wetting layer; ripening 
islands without wetting layer. 

equilibrium configuration is a uniform adsorbate layer on 
the substrate; this phase is known as the "Frank- Van der 
Merwe" phase. If it is less than zero, the equilibrium 
phase depends on the values of h min , L m i n , d m in, and 
Zmin which minimize AE. A few different situations can 
arise (we omit the "min" subscript to facilitate reading): 

1. If d — > oo, the islands undergo ripening, either di- 
rectly on the substrate (z — 0) or on a wetting layer 
(*^0); 

2. If d is finite and L < d, the equilibrium configura- 
tion is an array of islands either on the substrate 
( "Volmer- Weber" phase 22 ) or on a wetting layer 
( "Stranski-Krastanov" phased). 

3. If d is finite and L = d, the system is "cracked", 
that is, islands touch at their base. 

Following the example of Daruka et al^, we produce a 
set of phase diagrams in the a-8 plane. Figure 01 shows 
a typical phase diagram for specific values of the binding 
energies (eAA = 1 and esA = 1.1)- Such plots are col- 
lected in Fig. [S] for the whole array of values of €aa and 
esA investigated. 

We observe, first, that only one phase is possible when 
eAA > £sa (above the main diagonal of the plot), namely 
ripening islands without wetting layer. This is not sur- 
prising since it corresponds to the non-wetting case in 
which all terms of AE but ^interaction are negative. A 
more important conclusion that can be drawn from Fig. 
121 is that no set of control parameters leads to an array of 
islands as a stable configuration. This is in disagreement 
with the results of CJB 4 The discrepancy can be traced 
back to the choice of zq, the characteristic length for the 
decay of the adsorption energy at the surface, discussed 
in Section III B II The value of this parameter (0.39) was 
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FIG. 5: (Color online) Phase diagrams for all values of the 
binding energies considered; here ess = 1 and zq = 0.39. 
Each small image is a phase diagram in the a — 6 plane sim- 
ilar to that of Fig. 21 The colors correspond to the different 
phases: ^ Frank-van der Merwe; ^ ripening islands with 
wetting layer; f_H ripening islands without wetting layer. 

obtained, in the present model, by fitting to the LJ po- 
tential, while CJB chose it in an ad hoc manner. 



B. Extended model 
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FIG. 6: (Color online) Same as Fig. but for zq — 3.0. 
The various phases are as follows: Frank-van der Merwe; 
| | ripening islands with wetting layer; f_H ripening islands 
without wetting layer; ■ cracks; ■ Volmer- Weber. 

We may relax the constraint on the value of zq for 
a moment and set zq = 3, close to the value used by 
CJB, resulting in an increase of the influence of the sub- 
strate on the adatoms higher above the interface. The 
introduction of this new length scale in the problem al- 
lows some new features to appear. The corresponding 
phase diagrams are displayed in Fig. Consistent with 
CJB, we now observe the existence of two different stable 
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phases. In addition to the stable array of islands, we now 
find an equilibrium state consisting of islands touching at 
their base. Since this morphology is somewhat better de- 
scribed by a flat layer with very narrow troughs, we will 
refer to it as the "cracks" phase (shown in red in Fig.|BJ. 
The Stranski-Krastanov phase never appeared in the pa- 
rameter space we considered. 

It is not clear that potentials with zq = 3 do in fact 
exist, but it is certainly conceivable, in the case for in- 
stance of semiconductors, that the decay length of the 
surface energy is significantly larger than that of the LJ 
potential and, as a consequence, stable phases would ex- 
ist. Consideration of this extended system is consistent 
with our long term objectives since zq can be included as 
an adjustable parameter in the KMC calculations. 




Lattice mismatch (a, %) 



FIG. 7: (Color online) Phase diagram for the system with 
zo = 3.0, €aa ~ 1, and esA = 1.3. The various phases are: 
Frank- van der Merwe (FM), ripening islands with wetting 
layer (Rl), ripening islands without wetting layer (R2), cracks 
(C), and Volmer- Weber (VW). The dashed line indicates the 
smooth transition between the VW and the C phases. Inset: 
original phase diagram used to construct the main graph. 

In order to get more detailed information about the 
stable phases, we selected one of the phase diagrams of 
Fig. El (esA — 1-3, 6aa = 1) and repeated the calcula- 
tions on a finer grid (200x200). The resulting diagram 
is shown in Fig.0 Note that, in this graph, the layout of 
the phase boundaries have been drawn as a "guide to the 
eye" using the original data (shown in the inset); we do 
not know the analytical form (if it exists) of these curves. 
(The slightly "zigzagging" behaviour of the boundary is 
a consequence of allowing z to take only integer values 
in the minimization process.) We find that the transi- 
tion is not sharp between the Volmer- Weber (VW) and 
"cracks" (C) phases (hence the dashed line). Figurcdhas 
some features similar to the phase diagram computed by 
Daruka et al£, in particular the shape of the FM, Rl and 
R2 phases. 

The Frank-van der Merwe (FM) and the two ripening 
(Rl and R2) phases have no interesting intrinsic features; 
the first is flat, and the two others correspond to the 



minimizing parameters going to infinity. We are therefore 
more concerned about the two stable phases (C and VW). 
A characteristic quantity often measured is the island 
density n on the surface. In our model, since the center- 
to-center distance between islands is d, we simply have 
n = l I d. Figure |H1 shows how the density varies as a 
function of coverage and lattice misfit for the binding 
energies selected. In the C phase, this density is simply 
the cracks density, that is, the number of cracks per unit 
length. 

We have also computed the aspect ratio of the islands 
as a function of coverage and lattice mismatch; this is 
shown in Fig. The global behavior is as expected: 
as a increases, the elastic relaxation process becomes 
more and more efficient and the islands can afford an 
increase of their surface; this explanation holds for the 
cracks phase as well. 




Lattice mismatch (a, %) 

FIG. 8: (Color online) Island density 1/d in the two stable 
phases as a function of coverage 9 and lattice mismatch a. 
Here, zo = 3.0, eaa = 1, and esA = 1.3. 
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FIG. 9: (Color online) Aspect ratio h/L in the two stable 
phases as a function of coverage 9 and lattice mismatch a. 
Here zq = 3.0, €aa = 1, and esA = 1-3. 
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In both Figs. |Sland|51 the quantities shown by the color 
scale are conspicuously continuous at the VW-C bound- 
ary. This was to be expected since the very definition of 
these phases implies no jump in any quantity (VW phase 
for L < d, C phase for L = d). This boundary actually 
is the only second order phase transition; all others are 
first order. 

To our knowledge, a cracks phase has never been re- 
ported. We do not know if an equivalent 3D phenomenon 
has been observed experimentally. In our model, the oc- 
curence of such a phase is obviously related to the interac- 
tion energy term. Have we had assumed that interaction 
is infinite when the islands touch (as did CJB4), the C 
phase would have been entirely precluded. A stronger 
island-island repulsion might also lead to the appearance 
of a Stranski-Krastanov (SK) phase for intermediate cov- 
erage. While it remains to be verified, we expect this 
cracks phase to be very close, at finite temperature for 
instance, to a SK phase; this is however beyond the scope 
of the present work. 



IV. CONCLUSION 

We have presented a simple analytical model for the 
determination of the stable phases of strained heteroepi- 
taxial systems in (1 + 1) dimensions. The model was 
developed with a view of carrying out KMC simulations 
of the dynamics of the formation of islands. This is a 
very difficult task, but already we have made progress in 
this direction, on which we will report in a subsequent 
publication. In order for the present model to be ex- 
portable to a KMC code, all expressions were adjusted 
to an atomistic Lennard- Jones system. The present cal- 
culations reveal that, for parameters which arc consistent 
with the Lennard- Jones model, the array of islands is not 
a stable configuration of the system. If full consistency 
of the parameters is not imposed, and in particular if we 
relax the value of the decay length for the adsorption en- 
ergy (zo), then a stable array of islands arises. We argue 
that zq may in fact be viewed as an adjustable param- 
eter, which can be used to describe systems other than 
Lennard-Jones, in particular semiconductors. Our calcu- 
lations also reveal, in these conditions, the formation of 
a cracks phase — an array of islands touching at their 
base. While much work remains to be done, the present 
model is a first step in our aim of better understanding 
the formation of dislocation- free arrays of islands. 
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APPENDIX A: COMPUTATIONAL DETAILS 

We present here some details of the method used to 
fit the free parameters arising in the derivation of the 
continuous model presented in Section ITT1 

In what follows, the cutoff radius of the potential has 
been fixed to r c — 3.2, i.e., midway between the 5th and 
6th neighbour shells. The equilibrium distance between 
substrate atoms is set to 1 so that ass = 1/1.1119 = 
0.8993 (cf. Table UJ). In all calculations, the substrate 
has thickness between 50 and 100 layers, with the lower 
three maintained fixed to mimic the presence of the bulk. 
For every configuration, the energy was determined by 
relaxing the positions of the atoms using a conjugate- 
gradient algorithm. 



1. Surface Energy 

The parameters B and zq in Eq. i|13[l have been fitted 
to systems composed of 50 substrate layers and cover- 
age 8 £ {1, 2, 3, 4, 5, 10}. These configurations have been 
relaxed for a set of combinations of energy parameters 
(tAA,esA) £ {0.8, 0.9,1, 1.1, 1.2} 2 . Altogether, the min- 
imum total energies of 150 different systems have been 
computed. For every system, the difference between the 
total energy and the total energy of the equivalent system 
with esa = (-AA = ^SS has been calculated. The numeri- 
cal values for i5^ yer — E l g ycl in (|13fl have then been fitted 
to an equation of the form 

flayer _ flayer = ^(fl)^ _ ess) + C 2 (9)(e AA - e SA ). 

(Al) 

Figure ITUI shows the dependence of C\ and Ci on cover- 
age. C\ is completely determined (i.e., there are no free 
parameters), given by C\ = du$6, where uq = 3.364 is 
the cohesive energy for a cutoff radius of 3 (see Table [I}. 
C2 was fitted to a curve of the form B(l — e~ e / z °); the 
fit yields B = 2.53 and z = 0.39. 



2. Elastic Energy and Island-island Interaction 
Energy 

The parameters ai, 02, 61, 62, and c of Eqs. (|27|l and 
(|3l))l have been fitted to the same configurations as above. 
30 different islands have been generated, of width L in 
the range 20 to 80 (six values) and height h in the range 
1 to L (five values). Each of these configurations was 
relaxed with misfits of +1% and —1%. While the relaxed 
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FIG. 10: Comparison between the numerical values of Ci and 
Ci (in equation lAlfl and the theoretical curves (see text), (a) 
The curve has no free parameter, (b) The curve is the best 
fit to the data. 
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FIG. 12: Numerical data and theoretical curves for the strain 
energy given by J28I . 



displacements of the atoms of the substrate; for every sys- 
tem, the amplitude of the theoretical displacement field 
[which is close to, but not equal to Eq. (|19fl because of 
periodic boundary conditions] has been adjusted to the 
displacements of the substrate atoms. Figure ^] shows 
the agreement between the curves and the data. 

Figure^Jshows the good agreement between the strain 
energy obtained from the simulations and the analyti- 
cal expression Q2H[1. In Appendix iBl we elaborate on the 
choice of this equation. 



APPENDIX B: FUNCTIONAL FORM OF EQ. 
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FIG. 11: Numerical data and theoretical curves for K geo given 
by equation 1361 . 
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FIG. 13: Elastic energy for an island of given volume V as a 
function of its aspect ratio h/L. 



positions of the atoms depend on the sign of a, the energy 
does not, as could be expected. 

The numerical value of K geo has been found using the 



We mentioned in Section Hi B 21 that a better fit to the 
numerical data is obtained with a function of the form 
E(L,h) — Cv(L,h) [where v(L,h) has units of volume], 
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- Combe et al. 
Eq. (B3) 



W L 
Island height (h) 

FIG. 14: Elastic energy for an island of given width L as a 
function of its aspect height h. 



rather than E(L, h) = CVR(h/L). We show here a com- 
parison between our expression for v(L, H), Eq. (|2"7j) . and 
that used by CJB4; the latter is obtained by writing, first, 
the function R(h/L) in terms of the quantities used in the 
present paper: 



Rc.m(r)=A + Be- c ^-^\ 



with A = 0.13, B = 0.87 and C = -4.811 and r = h/L. 
Since v(L, h) = R(h/L)V, it is a simple matter to connect 
the two functional forms. We get 



V3 



h(L-h/2) (A + Be- Ch /( L - h M 



2 v ' ' V 
which must be compared with Eq. I|27|) . and 



R(r) 



1 



~cr I (1— r)\ 



cr(l - r/2) 



(B2) 



(B3) 



(Bl) 



Figure ED shows the the elastic energy of an island of 
fixed volume Fasa function of its aspect ratio, according 
to Eqs. I|B3|) and i|Bl|) . The two curves are different, but 
their general behaviour is very similar. Note in particular 
that the starting points coincide and that the end points 
are less than 2% apart. 

The situation is however very different for the elastic 
energy of an island of fixed width as a function of height 
(Fig. ITU) . CJB's expression for the energy yields a peak 
around h/L = 0.2 and has a minimum around h/L = 0.6; 
in contrast, our expression increases monotonically with 
height. Our preference for the form of Eq. (|28|) is mainly 
based on this observation and the quality of the fit to the 
LJ data (see Fig. 112)1 . 
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